Molecular subtypes, tumor microenvironment infiltration characterization and prognosis model based on cuproptosis in bladder cancer

Cuproptosis is a kind of cell death dependent on copper. We aimed to explore the functions of the cuproptosis in the tumor microenvironment (TME) and construct a cuproptosis-related prognosis signature in bladder cancer (BCa). Using BCa patients in the public cohort, the cuproptosis-related molecular subtypes and cuproptosis-related prognosis signature were developed. Three cuproptosis-related molecular subtypes, with different prognoses and TME characteristics, were identified in BCa. The cuproptosis-related prognosis signature can divide patients into high- and low-risk groups with different prognoses, TME characteristics, chemotherapeutic drug susceptibility and immunotherapeutic response. Low risk group patients had a favored prognosis and response to immunotherapy. The dysregulation of cuproptosis-related genes expression levels was validated in multiple BCa cells using in vitro experiments. Cuproptosis has an important role in the tumor progression and the characterization of TME in BCa. The cuproptosis-related prognosis signature is a useful biomarker that can reflect the prognosis, TME characteristics, immunotherapeutic response and chemotherapeutic drug susceptibility in BCa patients.


INTRODUCTION
Bladder cancer (BCa) is an aggressive cancer characterized by high rates of recurrence and metastasis (Babjuk et al., 2022;Witjes et al., 2021). The 2020 Global Cancer Statistics has reported that there were 573,278 new BCa patients worldwide, accounting for 3% of all newly diagnosed tumor patients; 212,536 deaths from BCa, accounting for 2.1% of all patients died of tumor (Sung et al., 2021) . Based on muscle-invasive status, BCa is broadly classified into nonmuscle invasive bladder cancer (NMIBC, 75%) and muscle invasive bladder cancer (MIBC, 25%) with distinct treatments and prognosis (Babjuk et al., 2022;Montironi et al., 2016;Witjes et al., 2021). Currently, targeted therapy, radiotherapy and immunotherapy for treating BCa patients, tumor recurrence and metastasis remain the main cause of treatment failure and fatal survival outcomes. Thus, identifying molecular subtypes of BCa is necessary to reduce the rates of recurrence and metastasis.
Based on gene expression, there are several types of molecular classifications, including TCGA (Robertson et al., 2018), Baylor (Mo et al., 2018), Lund (Lindgren et al., 2010), UNC (Damrauer et al., 2014), MDA (Choi et al., 2014), CIT (Rebouissou et al., 2014) and EUA (Kamoun et al., 2020).These molecular classifications can distinguish the biological characteristics and prognosis of BCa patients among molecular subtypes. However, the functions of these molecular classifications in treatment decision-making remain unidentified in BCa. Therefore, developing new molecular subtypes that can guide therapy need to be considered. Programmed cell death plays an important role in tumorigenesis, tumor progression and therapy (Labi & Erlacher, 2015;Lee et al., 2018;Li et al., 2019). Chemotherapeutic drugs that induce autophagy or apoptosis of cell death, like gemcitabine and cisplatin, have commonly been used in the treatment of BCa and achieved some clinical benefits (Goel et al., 2019;Lin et al., 2017). Cuproptosis occurs through the combination of copper with lipoylated components in the tricarboxylic acid cycle. This process promotes the loss of subsequentiron-sulfur cluster protein and the aggregation of the lipoylated protein, which leads to apoptosis (Tsvetkov et al., 2022). A tumor with a high concentration of lipoylated proteins is sensitive to cuproptosis, revealing that copper ionophore treatment may be an alternative treatment to these tumors that are resistant to traditional therapies.
We aimed to perform a comprehensive study to investigate the role of cuproptosis in BCa from a tumor microenvironment (TME) perspective. In this study, 1,306 BCa patients were classified into different cuproptosis-related subtypes, and the differences in prognosis and TME characteristics among the three subtypes were investigated. A cuproptosis prognosis signature was constructed to predict survival outcomes and characterize the chemotherapeutic drug susceptibility and the immune landscape of BCa. The results revealed that the cuproptosis score may be a useful prognostic marker. In addition, we further performed in vitro experiments to explore the expression pattern of cuproptosis genes and validate the performance of the cuproptosis-related prognosis signature. Figure S1 demonstrated the process of this study.

Data collection
In Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), datasets were included based on the following criteria: (1) Histologically diagnosed as BCa; (2) Availability of RNA expression data; (3) Availability of prognostic information; (4) More than 50 tumor samples. In this way, the normalized gene expression data of five BCa datasets (GSE48277, GSE32894, GSE31684, GSE48075 and GSE13507) with relevant prognostic data were obtained from the GEO. The RNA-sequence data (fragments per kilobase million, FPKM) with clinical information of patients were obtained from a BCa dataset in The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov). The gene expression data and corresponding clinicopathological and follow-up information of metastatic urothelial cancer patients who accepted atezolizumab (PD-L1 inhibitor) were obtained from the IMvigor210 trial (http://research-pub.gene.com/IMvigor210CoreBiologies). BCa patients in five GEO datasets and IMvigor210 trial were combined through the ''Combat'' method from the ''sva'' R package (https://doi.org/doi:10.18129/B9.bioc.sva) (Leek et al., 2012) to minimize the batch effects. ''Combat'' method could adjust data for batch effects through an empirical Bayes framework and return an expression matrix that has been corrected for batch effects. As a result, there were 1,306 BCa patients in the combined dataset Table S1 showed the main information on the combined dataset.

Associations of molecular subtypes with clinicopathological features
The correlation of molecular subtypes with clinical features and prognoses was investigated. In addition, Kaplan-Meier and log-rank tests were performed using R packages ''survminer'' and ''survival''.

Cuproptosis gene cluster and cuproptosis-related prognosis signature
We further tried to construct the cuproptosis gene cluster and a cuproptosis-related prognosis signature. First, the univariate Cox regression and Kaplan-Meier analysis were performed on DEGs to obtain the prognostic DEGs. Second, Patients were divided into distinct cuproptosis gene clusters according to the gene expression of the prognostic DEGs using the R package ''ConsensusClusterPlus''. Finally prognostic DEGs were introduced into the least absolute shrinkage and selection operator (LASSO) model was used to develop the cuproptosis-related prognosis signature via the R package ''glmnet'', and the penalty regularization parameter lambda (λ) was chosen using 10-fold cross validation to obtain genes with non-zero coefficients and minimize the mean square error. Meanwhile, the minimal λ was identified to obtain the genes. The risk scores of the cuproptosis-related prognosis signature were calculated as follows: cuproptosis-related prognosis signature = Coef 1 × expression of gene 1 + Coef 2 × expression of gene 2 + . . . + Coef m × expression of gene m. Coef denotes the corresponding coefficient of the gene.

Chemotherapeutic drugs susceptibility analysis and mutation
The semi-inhibitory concentration (IC50) of BCa patients was evaluated based on the R package ''pRRophetic''. The somatic mutations of BCa patients in TCGA were demonstrated by the R package ''maftools''.

Development of the nomogram
The cuproptosis-related prognosis signature and clinical factors were subjected to univariate cox analysis. Factors significant in univariate cox analysis were then subjected to multivariate cox analysis to develop a nomogram. Decision curve analysis (DCA), Time-dependent receiver operating characteristic curve (ROC) curves and Calibration plots were applied to explore the performance of the nomogram and clinical factors.

Cell culture and Real-time quantitative PCR (RT-qPCR)
All cell experiments were approved by the Ethics Committee of Shanghai Tenth People's Hospital (approval number 2021KN108).
BCa cell lines (5637, T24, RT4, J82, and UMUC3) and human normal bladder epithelial cell line SV-HUC-1 were all purchased from the Chinese Academy of Sciences. RT4 was cultured in McCoy's 5A medium (Thermo Fisher Scientific, Waltham, MA, USA); SV-HUC-1 was cultured in F12k medium (Sigma-Aldrich, St. Louis, MO, USA); J82 was cultured in DMEM (Thermo Fisher Scientific); T24, 5637 and UMUC3 were cultured in RPMI-1640 medium (Thermo Fisher Scientific). All cells were cultured at 37 • C and 5% CO2. All cell culture medias were supplemented with 10% fetal bovine serum (Gibco, Waltham, MA, USA). The protocols of total RNA extraction and RT-qPCR have been described previously (Wang et al., 2021). Data were analyzed with Bio-Rad software (Bio-Rad, Hercules, California, USA). Data were normalized by β-actin expression level in each sample. 2 − Ct method was used to calculate relative expression changes (Livak & Schmittgen, 2001). Primer sequences are presented in Table S2.

Statistics
Principal component analysis (PCA) was conducted using the R package ''pcaMethods''. The correlations among factors were investigated using the R package ''corrplot'', ''circlize'', ''igraph'' and ''reshape2''. The cutoff points of the cuproptosis score and risk score were determined by the R package ''survminer''. SPSS 23.0 (IBM, NY, USA) and R v.4.0.3 software (R Core Team, 2020) and were used for statistical analysis. Variables between subtypes or groups were compared using Mann-Whitney U test, Chi-square test, or Student's t -test, as appropriate. A P < 0.05 in two-tailed analyses was statistical significance, with a Bonferroni correction applied for pairwise comparisons. Figure 1A showed the chromosome locations of the cuproptosis-related genes. GeneMANIA network built around cuproptosis-related genes demonstrated the top 20 genes interacted with cuproptosis-related genes ( Fig. 1B). High interactions within cuproptosis-related genes were also showed in this network. The correlation network among cuproptosis-related genes were shown in Fig. 1C. Then, we investigated somatic mutations in cuproptosis-related genes in TCGA and discovered relatively low somatic mutations in these genes (Fig. 1D). Among them, CDKN2A had the highest somatic mutations while others had five or less than five somatic mutations.

Identification of cuproptosis subtypes in BCa
Through the ''Combat'' algorithm, the batch effects were eliminated in the combined dataset (including GSE13507, GSE31684, GSE32894, GSE48075, GSE48277 and IMvigor210 trial) (Figs. S2A, S2B). The gene expression data of cuproptosis-related genes were used to categorize the BCa patients in the combined dataset. As a result, K = 3 was an optimal value to classify the BCa patients in the combined dataset into Cluster 1/2/3 (n = 515/427/364) (Figs. S2C-S2E). PCA plot showed novel differences among the three clusters (Fig. 1E). Survival analysis revealed a poorer overall survival (OS) in Cluster 2/3 patients than those in Cluster 1 patients (P = 0.027; Fig. 1F, Fig. S2F). Ten cuproptosis-related genes had significantly different expression patterns among these three clusters (Fig. 1G, Fig. S2G). In addition, Cluster 2 and Cluster 3 were preferentially related to the high T stage and higher risk of death (Fig. 1H). (G) Differences in clinicopathologic features between three distinct subtypes. BCa, bladder cancer; OS, overall survival; NMIBC, nonmuscle invasive bladder cancer; MIBC, muscle invasive bladder cancer; TCGA, the Cancer Genome Atlas. *, p < 0.01, * * * , p < 0.001. p value was calculated with t test, except survival analysis was analyzed using a two-sided log-rank test.

Characteristics of cuproptosis score and TIICs infiltration in the cuproptosis subtypes
Survival analysis showed that cuproptosis score was a protective biomarker in BCa patients ( Fig. 2A). We then explored the correlation of the cuproptosis score with ten cuproptosisrelated genes and cuproptosis subtypes. The results revealed that ten cuproptosis-related genes all had a positive correlation with the cuproptosis score (Fig. 2B, Fig. S3). Cluster 1 patients had the highest cuproptosis scores while Cluster 2 patients had the lowest cuproptosis scores (Fig. 2C).
We further investigated the differences in the proportion of TIICs among cuproptosis subtypes. The relative proportion of the 22 TIICs in each BCa patient was calculated and showed in Figs. 2D, 2E. Eighteen of 22 TIICs had significantly differently composition among cuproptosis subtypes (Fig. 2E), which revealed the diverse tumor immune cell microenvironments among cuproptosis subtypes.

Identification of DEGs and cuproptosis gene cluster
As patients in Cluster 2/3 had similar survival outcomes, we combined Cluster 2 and Cluster 3 (named Cluster 2/3) and identified cuproptosis subtype-related DEGs between Cluster 2/3 and Cluster 1. In this way, 291 cuproptosis subtype-related DEGs were obtained (Fig. 3A), and DEGs that were up-regulated in Cluster 2/3 were enriched in many cancer-promoting biological processes (Figs. 3B, 3C), including Cell division, Cell cycle, DNA biosynthetic process, DNA replication, DNA replication initiation. The cuproptosis subtype-related DEGs were then introduced to univariate Cox analysis and Kaplan-Meier analysis in TCGA and combined dataset to obtain the prognostic DEGs. In this way, 71 prognostic DEGs were the overlapping genes related to OS in TCGA and combined data (Fig. 3D, Table S3). Based on the above genes, BCa patients in the combined dataset were divided into two genomic subtypes using the consensus clustering algorithm, namely GeneCluster 1 and GeneCluster 2 (Figs. S4A-S4F). PCA plot showed significant differences in the cuproptosis transcription profiles among two cuproptosis gene subtypes (Fig. 3E). Patients in GeneCluster 2 had significantly poorer OS than those in GeneCluster 1 (P < 0.001, Fig. 3F), and BCa patients in GeneCluster 2 had significantly higher cuproptosis scores compared with BCa patients in GeneCluster 1 (Fig. 3G). In addition, there were significant differences in clinicopathological features (Fig. S4G), TME characteristics (Figs. S5A, S5B) and cuproptosis-related gene expression (Fig. 3H) between GeneCluster 1 and GeneCluster 2.

Development of the cuproptosis-related prognosis signature
Based on the minimal λ (−2.840), 16 hub genes with their corresponding coefficients were obtained among prognostic DEGs to develop the cuproptosis-related prognosis signature (Fig. 4A, Fig. S6A), in which the partial likelihood deviance was at its minimum (Fig.  4B). The correlations of 16 hub genes in the cuproptosis-related prognosis signature were demonstrated in Fig. S6B, suggesting that these genes were not highly correlated with each other. The risk scores were calculated by summing the genes weighted by corresponding coefficients, which was mentioned above. The optimal cutoff point for risk scores was 0.812 determined by the R package ''survminer'' in TCGA. BCa patients were classified into high-and low-risk groups according to the optimal cutoff point, and high-risk patients had significantly lower survival rates and poorer OS than low-risk patients (Figs. 4C, 4D). Based on the prognostic DEGs, PCA was performed and revealed that BCa patients in low or high risk groups were distributed in discrete directions indicating novel differences in the prognostic DEGs between the high-and low-risk patients (Fig. 4E). To validate the prognostic value of the cuproptosis-related prognosis signature, BCa patients in the combined dataset were also divided into low-and high-risk groups. Kaplan-Meier curve showed that patients in high risk group also had poor OS than those in the low risk group in the combined dataset (Fig. 4F). The alluvial diagram illustrated the distribution of BCa patients in three cuproptosis subtypes, two risk groups, two gene subtypes, two cuproptosis score groups and survival outcomes (Fig. 4G). Patients in Cluster 2/3 and GeneCluster 1 had higher risk scores compared with those in Cluster 1 and GeneCluster 2, respectively (Figs. 4H, 4I). In TCGA, high-risk patients were preferentially correlated with higher age, higher TNM stage and low survival rate (Fig. 4J, Figs. S6C-S6G).

Characteristics of the signature
We compared somatic mutations, TMB, cuproptosis score and IC50 values of chemotherapeutic drugs between low-and high-risk patients. In TCGA, the top ten mutated genes were demonstrated in Figs. 5A and 5B, respectively. High-risk patients had lower TMB than those in low risk patients (P = 0.03, Fig. 5C). There was a negative correlation between TMB and risk scores (Fig. 5D, R =−0.17, P < 0.001). Similarly, in the combined dataset, low risk patients had higher cuproptosis score compared with high-risk patients (P < 0.001, Fig. 5E), and risk score was negatively related to the cuproptosis score (Fig. 5F). In addition, among seven common chemotherapeutic drugs used in BCa, the IC50 of five drugs were significantly different between two risk groups (Fig. 5G). GSEA revealed that several cancer-promoting pathways were highly enriched in high risk patient, including 'Focal adhesion', 'MAPK signaling pathway' and 'Pathway in cancer', while the 'Positive regulation of necrotic cell death' was enriched in low risk BCa patients (Fig. 5H). GSVA showed similar results with the 'MAPK signaling pathway', 'JAK STAT signaling pathway', 'Focal adhesion' and 'Pathway in cancer' and 'WNT signaling pathway' enriched in high risk BCa patients (Fig. 5I).

Characteristics of the TME in distinct risk groups
Fourteen of 22 TIICs had different compositions between two risk groups (Fig. 6A, Table  S4), which revealed the diverse tumor immune cell microenvironments between risk groups. Specifically, the proportions of macrophages were higher in high-risk patients compared with low risk patients, while T cells follicular helper, T cells CD8, B cells memory and monocytes contained significantly lower infiltration levels in high-risk patients compared to those in low risk patients. In addition, the risk score had a significant correlation with 12 TIICs (Fig. 6B, Fig. S7). We then investigated the relationship between the expression of some immune checkpoints and the risk score. The results revealed that some immune checkpoints had different expression levels between two risk groups, including CD274, PDCD1 and CTLA4 (Fig. 6C), and risk score had a correlation with some immune checkpoints (Fig. 6D). High-risk patients had higher immune scores, stromal scores and TIDE scores compared with low-risk patients (Figs. 6E, 6F). Based on the TIDE scores, high-risk patients had a lower proportion of responders compared to those in low-risk patients if received immunotherapy (Figs. 6G, 6H). These results indicated that this model could be a useful predictor of the efficacy of immunotherapy. GSEA and GSVA revealed that high-risk patients had multiple immune-related pathways enriched, TCGA, the Cancer Genome Atlas; OS, overall survival. p value was calculated with t test, Bonferroni correction applied for pairwise comparisons. Survival analysis was analyzed using a two-sided log-rank test.

Construction of the nomogram
To overcome the inconvenient clinical utility of the signature (risk score), a cuproptosisrelated prognosis signature-based nomogram was constructed to predict the survival outcomes of BCa patients. The results of univariate and multivariate Cox regression analysis showed that risk score and T stage remained the independent prognostic factors (Fig. 7A), so the signature and T stage were selected to develop the nomogram (Fig. 7B).
The histogram also showed that patients with T stage ≥T3 and high-risk patients had higher death rates than others (P < 0.001, Fig. 7C). The nomogram achieved the optimal performance in the prediction of patients' survival outcomes compared with other factors (Figs. 7D-7H). Furthermore, the nomogram achieved higher clinical net benefit than the T stage and risk group (Figs. 7I-7K).

The experimental validation of the cuproptosis-related genes
The mRNA expression pattern of the ten cuproptosis-related genes in five human BCa cell lines (T24, 5637, RT4, J82, and UMUC3) and human normal bladder epithelial cell line SV-HUC-1 was determined using RT-qPCR. The results showed significant differences in the expression levels of nine cuproptosis-related genes (except LIAS) between BCa and normal bladder epithelial cells (Fig. 8).

DISCUSSION
The heterogeneity of BCa leads to highly variable survival outcomes and makes antitumor treatment a demanding task. It is reported that the response rates to neoadjuvant chemotherapy were less than 50% in BCa (Grossman et al., 2003;Kitamura et al., 2014;Zargar et al., 2015), and the response rates to immunotherapy were only about 25% in BCa. (Bellmunt et al., 2017;Sharma et al., 2017;van der Heijden et al., 2021). Thus, it is necessary to accurately identify the subtypes of BCa for personalized therapy. Previous studies have developed various molecular classifications of BCa (Kamoun et al., 2020;Robertson et al., 2018;Sjödahl et al., 2012), but the molecular subtypes of these molecular classifications still cannot be used to determine the treatment strategy due to the considerable heterogeneity among the subtypes. Therefore, a more accurate molecular subtype of BCa is needed to guide personalized therapy and improve survival outcomes. Cuproptosis is a newly discovered cell death, and it depends on the combination of copper with lipoylated TCA cycle proteins (Tsvetkov et al., 2022). As copper is one of the cofactors for essential enzymes in human cells (Kim, Nevitt & Thiele, 2008), the identification of molecular subtypes based on cuproptosis may be a useful direction for clinical treatment.
To understand the significant role of cuproptosis in BCa, our study established three subtypes in BCa patients based on ten cuproptosis-related genes. Significant differences in prognoses and biological characteristics were found among the three subtypes. Several cancer-promoting biological pathways were highly enriched in Cluster 2/3, which partially explained the poor survival outcomes in Cluster 2/3. The proportions of the TIICs also differed significantly among the three subtypes. To investigate the functions of cuproptosis in BCa patients' prognosis, we used the GSVA algorithm to create a biomarker named cuproptosis score according to ten cuproptosis-related genes. The results revealed that a high cuproptosis score was correlated with a favored prognosis, and BCa patients in Cluster 1 characterized by a favored prognosis have the highest cuproptosis score among the three subtypes. In addition, two gene subtypes (GeneCluster 1 and GeneCluster 2) were obtained based on the prognostic DEGs among the three subtypes. Significant differences in the RNA expression levels of most cuproptosis-related genes were found between GeneCluster 1 and GeneCluster 2, and BCa patients in GeneCluster 1 have significantly lower cuproptosis scores and poor OS. These findings indicated that cuproptosis has an important role in the tumor progression, TME and prognosis in BCa, and BCa patients with different cuproptosis patterns had different TME characteristics and prognoses.
Many studies developed prognosis signatures to predict the prognoses of tumor patients based on various kinds of cell death, including pyroptosis and ferroptosis (Song et al., 2021;Sun et al., 2021;Zhang et al., 2021). However, most of these prognosis signatures are based on single dataset without external validation or experimental validation, which limits the reproducibility and generalizability of the prognosis signatures. We aimed to established a useful prognosis signature and validated its prognostic value in other datasets. In addition, experimental validation was performed which could be lacked in other studies. In this study, a cuproptosis-related prognosis signature based on 16 prognostic genes was developed. As predicted, the cuproptosis-related prognosis signature was related to patients' prognoses, high risk score was related to advanced TNM stage and poor survival outcomes. BCa patients in Cluster 1 and GeneCluster 2 characterized by favored prognosis have significantly lower risk scores compared to Cluster 2/3 and GeneCluster 1, respectively. The prognostic value of the cuproptosis-related prognosis signature was further explored in the combined dataset for external validation. As predicted, the cuproptosis-related prognosis signature still had a novel performance in survival prediction in the combined dataset, indicating that the cuproptosis-related prognosis signature may be a novel prognostic biomarker. Additionally, patients in low risk group have significantly higher IC50 values of five chemotherapeutic drugs compared to those in high-risk patients, which revealed that the cuproptosis-related prognosis signature may predict the therapeutic response of chemotherapeutic drugs and high-risk patients may benefit from chemotherapy. GSEA and GSVA further showed that there were several cancer-promoting pathways highly enriched in high-risk patients and a pathway of 'positive regulation of necrotic cell death' enriched in low-risk patients, revealing that two risk groups divided by the cuproptosis-related prognosis signature have distinct biological characteristics, thereby contributed to differences in prognosis and therapeutic response between two risk groups. The cuproptosis-related prognosis signature may be useful in survival prediction and personalized therapy in BCa.
While the use of chemotherapeutic drugs cannot produce a durable response and may lead to some serious side effects, immune checkpoint inhibitors (ICIs) is playing a vital role in the therapy of BCa patients. Clinical research has proved the safety of ICIs in human tumors (Rotte , 2019). For advanced and metastatic urothelial carcinoma in the post-platinum setting, several ICIs have been suggested as the second-line treatment (Bellmunt et al., 2017;Patel et al., 2018;Powles et al., 2017;Rosenberg et al., 2016;Sharma et al., 2017;van der Heijden et al., 2021). However, only about 25% BCa patients were sensitive to immunotherapy (Bellmunt et al., 2017;Sharma et al., 2017;van der Heijden et al., 2021), which highlight the necessary of constructing a reliable signature to identify subgroups of BCa patients benefit from ICIs. However, the specific mechanism between cuproptosis and TME is not clear. We then explored the function of cuproptosis in TME. Based on the cuproptosis-related prognosis signature, two risk groups had novel differences in the proportions of 14 TIICs and expression patterns of 27 immune checkpoints, and risk scores were highly related to most of the TIICs and immune checkpoints. Specifically, high-risk patients have a significantly higher abundance CD8+ T cells, B cells memory, CD4+ naive T cells and T-follicular helper cells, and low-risk patients have a significantly higher abundance of M2 macrophages and neutrophils. CD8+ T cells play an important role in tumor cell killing and therapeutic sensitivity of ICIs (Fridman et al., 2017;Powles et al., 2014), and patients with high abundance of tumor-infiltrating T cells have a favored prognosis in BCa (Sharma et al., 2007). Previous studies also showed that B cells in TME were correlated with better prognosis in BCa (Zirakzadeh et al., 2020), colorectal cancer (Berntsson et al., 2016) and sarcoma (Petitprez et al., 2020). In this study, T cells and B cells were highly enriched in the low-risk patients, which not only explained the good prognosis of low-risk patients but also revealed that low-risk patients possibly may benefit from immunotherapy. On the contrary, M2 macrophages can induce tumor immunosuppression and promote tumor progression through releasing proinflammatory cytokines including IL-10/12/13/16, INF-γ and transforming growth factor-β (Chen et al., 2019;Jeannin et al., 2018;Sica et al., 2006;Yunna et al., 2020), and tumor-associated neutrophil was another important member of TIICs and can release several factors leading to tumor progression (Eruslanov, 2017;Galdiero et al., 2013;Mantovani et al., 2011). Consistent with previous studies, an increased abundance of M2 macrophages and neutrophils was found in high risk group. High-risk patients have significantly higher TME scores than those in low-risk patients, which revealed that low-risk patients have lower relative infiltration of stromal cells and immunocytes in TME. As for the TIDE score, low-risk patients have higher responder rates and lower TIDE scores. TIDE was regarded as a surrogate biomarker that can predict response to ICIs, and a lower TIDE score is related to better ICIs response and better prognosis (Jiang et al., 2018). In our study, low TIDE scores contribute to the relatively high ICIs response and good prognosis in low-risk patients. GSEA further validated the dysfunction of the immune system in high-risk patients. These results illustrate that there were distinct immunophenotypes between low-and high-risk groups, and the cuproptosisrelated prognosis signature is associated with the immunotherapy response. In this way, the cuproptosis-related prognosis signature may have the potential to facilitate individualized immunotherapy in BCa. As the cuproptosis-related prognosis signature could divided BCa patients into subgroups with different prognosis and immunophenotypes, we developed a nomogram combining the cuproptosis-related prognosis signature with the T stage to facilitate the integration of the cuproptosis-related prognosis signature into clinical practice. The ROC curves, Calibration plots and DCA revealed that the nomogram was a useful evaluation tool for prognostic predictions. This nomogram facilitated the use of the cuproptosis-related prognosis signature and improved the performance of survival prediction.
We further explore the mRNA expression pattern of the ten cuproptosis-related genes in human BCa cell lines using qPCR. The results revealed that most of cuproptosis-related genes exhibited dysregulated expression between a human normal bladder epithelial cell line and BCa cell lines, suggesting that dysregulated expression levels of cuproptosis-related genes served an important role in the carcinogenesis of BCa. Specifically, CDKN2A was downregulated in four BCa cell lines compared with the human normal bladder epithelial cell line. CDKN2A (also called p16) was an important tumor suppressor that prevents carcinogenesis through induction of senescence cell growth and arrest (Rayess, Wang & Srivatsan, 2012). CDKN2A promoter methylation is a frequent epigenetic event, which could result in silencing of CDKN2A expression and thereby leading to uncontrolled cell proliferation and carcinogenesis in several tumors (Riese et al., 1999;Samowitz et al., 2005;Saridaki et al., 2003) including BCa (Worst et al., 2018). However, the tumor suppression mechanism of CDKN2A from a cuproptosis perspective had not been investigated, which is worthy of studying.
This study has some limitations. First, additional experimental verification is needed to validate the relationship between TME and cuproptosis. Furthermore, some clinical factors including therapy were unavailable, which may lead to different cuprotosis state and TME.
In conclusion, we identified three subtypes of cuproptosis-associated molecules in BCa with distinct prognoses. Cuproptosis-related genes play important role in TME heterogeneity in BCa. The cuproptosis-related prognosis signature can distinguish the survival outcomes, biological characteristics, TIICs infiltration and immunotherapeutic response in BCa patients.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
The authors received no funding for this work.